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Role of isospin in the nuclear liquid-gas phase transition 
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We study the thermodynamics of asymmetric nuclear matter using a mean-field approximation 
with a Skyrme effective interaction, in order to establish its phase diagram and more particularly 
the influence of isospin on the order of the transition. A new statistical method is introduced to 
study the thermodynamics of a multifluid system, keeping only one density fixed, the others being 
replaced by their intensive conjugated variables. In this ensemble, phase coexistence reduces to 
a simple one-dimensional Maxwell construction. For a fixed temperature under a critical value, a 
coexistence line is obtained in the plane of neutron and proton chemical potentials. Along this line, 
the grand potential presents a discontinuous slope, showing that the transition is first order except 
at the two ending points where it becomes second order. This result is not in contradiction with 
the already reported occurrence of a continuous transformation when a constant proton fraction is 
imposed. Indeed, the proton fraction being an order parameter in asymmetric matter, the constraint 
can only be fulfilled by gradual phase mixing along the first-order phase-transition line, leading to 
a continuous pressure. 

PACS numbers: 24.10.Pa,64.60.Fr,68.35.Rh 


I. INTRODUCTION 

Nucleons in atomic nuclei interact through a finite- 
range attractive and a short-range repulsive force. For 
systems of particles interacting this way, one expects to 
find a phase transition analogous to the liquid-gas tran¬ 
sition of a van der Waals fluid 0. As a matter of fact, it 
is recognized that symmetric nuclear matter should un¬ 
dergo a first-order transition between a low (gas) and a 
high (liquid) density phase up to a critical temperature 
|2; sj 13 ■ For such a system containing an equal number 
of neutrons and protons, isospin symmetry imposes that 
the nucleons behave as one single fluid and one expects a 
discontinuous density versus pressure equation of state. 

The case of asymmetric matter is more complex to 
study, since there is an additional degree of freedom 
to consider: the isospin. Such matter plays an impor¬ 
tant role in astrophysics, where neutron-rich systems are 
involved in neutron stars and type-II supernova evolu¬ 
tion For asymmetric systems containing a fixed 

proton fraction, it has been shown that the thermody¬ 
namic transformations result in a continuous evolution 
of the observables. In particular, the system density is 
a continuous function of the pressure. This has been 
interpreted as the occurrence of a continuous transition 
Hi i in mm. We will show in the present article 
that this conclusion is not correct: it results from a con¬ 
fusion between the notions of ’continuous transition’ and 
’continuous transformation’. Indeed, the phenomenon 
of isospin distillation demonstrates that the proton frac¬ 
tion is an order parameter in asymmetric nuclear matter. 
Thus, when the proton fraction is kept constant, the sys¬ 
tem is forced to follow the first-order phase-transition 
line, hiding the discontinuity of the thermodynamic- 
potential first derivative. 

The plan of the paper is as follows: we calculate in 
Section II the nuclear matter grand potential in the 


mean-field approximation using a Skyrme (Sly230a) en¬ 
ergy density functional (T3 | . The thermodynamics of 
nuclear matter has been addressed earlier with differ¬ 
ent effective forces El 13 m m The construction 
of phase coexistence is presented in section III. It con¬ 
sists in correcting the mean-field instabilities with the 
introduction of phase separation, in order to build the 
concave envelope of the thermodynamic potential. To 
perform this Gibbs construction with an arbitrary num¬ 
ber of conserved quantities, we introduce a new method 
that reduces this multidimensional problem to a simple 
one-dimensional Maxwell construction on a carefully de¬ 
fined statistical potential. The results of this analysis 
are presented in section IV. The grand-canonical poten¬ 
tial presents a discontinuous derivative on both sides of 
a bi-dimensional manifold limited by a critical line in 
the 3-dimensional space including the temperature and 
the proton and neutron chemical potentials. The transi¬ 
tion is thus first order for all associated proton fractions. 
Only the critical line corresponds to continuous transi¬ 
tions. An interesting point is that in this 3-dimensional 
problem the critical line can be characterized by addi¬ 
tional critical exponents as the chemical potential ap¬ 
proaches its critical value for a given temperature. We 
first present the coexistence region obtained for a fixed 
temperature, stressing the isospin properties of phase co¬ 
existence. Then, the thermodynamic consequences of a 
transformation at constant proton fraction are analyzed. 
We show that this transformation forces the system to 
follow the coexistence line, and this is the generic behav¬ 
ior expected when a conservation law acts on an order 
parameter. In this case, the first-order phase transition 
results in a continuous transformation from a diluted to a 
dense system through a phase coexistence, which should 
not be confused with a continuous transition. Indeed, 
one should clearly make the difference between a trans¬ 
formation, which is a specific path in the space of ther- 
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modynamic variables, and a phase transition which is an 
anomaly of the thermodynamic potential considered in 
the total space of thermodynamic variables . The 
temperature dependence of the phase diagram is finally 
explored, and the critical behaviour is presented. 


II. THERMODYNAMICS WITH SKYRME 
FORCES 

A. Isospin dependent energy functional 

We study the case of nuclear matter in a mean-field 
approach, with the SLy230a Skyrme effective interaction 
) 1 . j | . This local interaction allows to introduce an energy 
density H(r) so that the total energy for a system of 
nucleons in a Slater determinant | xf> > reads : 

(V'l = j H(v)dv (1) 

where H is the Hamiltonian of the system. The energy 
density H is a functional of the particle densities p q and 
kinetic densities r q for neutrons (q = n) and protons 
(q = p). Denoting p q the one-body density matrix of 
the particles of type q , those quantities are expressed as 
follows: p q {r) = (r|p,|r) and r q (r) = (r\^pp q p\r), such 

j-2 

that —T q is the kinetic energy density. 

For later convenience it is useful to introduce the 
isoscalar and isovector densities : 

P = Pn + Pp > T = T n + Tp 
P3 = Pn Pp ? 73 = 7"n 7~p 

In the case of homogeneous, spin-saturated matter 
with no coulomb interaction, four terms contribute to 
the energy density : 

H = JC + Ho+H3 +H ef f (3) 

In this expression, K, is the kinetic-energy term, Ho 
a density-independent two-body term, H 3 a density- 
dependent term, and H e f / a momentum-dependent term: 



h 2 

(4) 

/c 

2 TO 

Ho 

= Cop~ + D 0 P 3 

(5) 

n 3 

= C 3 p °+ 2 + d 3P ° pi 

(6) 

Heff 

= CeffPT + D effP 3 T 3 

(7) 


The coefficients Ci and Di , associated respectively 
with the symmetry and asymmetry contributions, are lin¬ 
ear combinations of the traditional Skyrme parameters : 

Co = 3fo/8 

Do = — to(2xo 1)/8 

C 3 = t 3 /16 / o', 

D 3 =-t 3 ( 2^ 3 + l)/48 

C e ff = [3fi + < 2 ( 4^2 + 5)]/16 

D e ff = [< 2 ( 2 x 2 + 1) — fi(2xi + 1)]/16 


To illustrate the SLy energy functional we present in 
figme^the energy density and the energy per particle as 
a function of the total particle density for various proton 
fractions. The kinetic-energy term has been computed 
integrating over the Fermi spheres associated with the 
considered proton and neutron densities. The minimum 
of the energy per particle is the saturation point of sym¬ 
metric matter po = 0.16 fm~ 3 and Eq = —15.99 MeV, 
while the curvature gives the incompressibility I\ — 230.9 
MeV. The Z/A = 0 curves correspond to pure neutron 
matter which does not saturate, the Sly forces b eing fit¬ 
ted on realistic neutron-matter EOS calculations fl3 . [l8j . 



FIG. 1: Sly230a functional for infinite nuclear matter : en¬ 
ergy per particle (upper part) and energy density (lower part) 
as functions of the total nucleon density for regularly spaced 
proton fractions from Z/A=0.5 to Z/A=0, in steps of 0.05. 

The mean-field effective Hamiltonian W q for each 
particle-type q is defined by the relation 8(H) = 

Tr(W q 5p q ) i.e. as a functional derivative of the energy 
density. In the case of Skyrme interactions, this leads to 
the expression : 


W q =™*- 
q dr q h 2 


an 

9Pq 


2 m: 


:f + U q 


(9) 


where U q = d Pq H is the local mean-field potential : 


u q = U 0q + U 3q + Ueff q (10) 

= [2Cop + (u + 2)C 3 p a+1 + a D 3 p a 1 p\ + C e //r] 
± [2Dop 3 + 2D 3 p a p 3 + D eff T 3 \ 


and the effective mass m* is defined by : 


h 2 

2m* 


2 m 


+ CeffP ± D eff p 3 


( 11 ) 
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In both expressions, the ± sign refers to neutrons (+) or 
protons (—). 

In uniform nuclear matter, the eigenstates of this 
mean-field Hamiltonian are spin-up or spin-down plane 
waves with momenta p;. The single-particle energies are 
given by : 




+ U q 


( 12 ) 


Equations d and d define a self-consistent prob¬ 
lem since m* depends on the densities according to 
eq. EJ. For each pair (/ Z n ,p ' p ) a unique solution ( p n ,p p ) 
is found by iteratively solving the self-consistency be¬ 
tween p rlJ) and m^ p . Then eq. m is used to calculate 
T n , P . These quantities allow to compute the one-body 
partition sum : 

Z 0 = Tr[e-P^ n +w p -p n N n -p, P N p )} = Z u Z p ( 2 0) 


B. Finite temperature 


where each partition sum Zq can be expressed as a func¬ 
tion of the corresponding kinetic energy density : 


Within a mean-field approach, thermodynamic rela¬ 
tions are easy to derive in the grand-canonical ensem¬ 
ble. For a system of neutrons and protons, the grand- 
canonical constraint imposes the average value of three 
observables: energy, number of protons and neutrons. 
Maximizing the Shannon entropy with these three con¬ 
straints leads to an equilibrium partition sum : 

Z GC = Tr[e-^ +anfln+a ^ ) ] (13) 

where the inverse temperature [3 = 1/kT is the La¬ 
grange multiplier associated with the energy constraint, 
and a q = (3p q are the Lagrange multipliers controlling 
the particle numbers < N q >, p q being the chemical po¬ 
tentials. The Lagrange parameters fulfill the equations 
of state : 


(H) = -fyln Z GC (14) 

(N q ) = d aq \n Z GC (15) 


The self-consistent mean-field approximation amounts 
to use independent particle states as trial density ma¬ 
trices in the maximum-entropy variational principle: the 
single-particle states of energy e q are then occupied ac¬ 
cording to the Fermi-Dirac distribution |ij| : 


q 1 + exp(/?(e* - p q )) 


(16) 


In infinite matter, n* is a continuous distribution n q (p) 
and the densities p q and r q read : 

r°° Anp 2 

Pq = 2 J n q (p)-^-dp (17) 

f°° 7} 2 Attv 2 

Tq = 2 J 0 ( 18 ) 

where the factor 2 come from the spin degeneracy. 

The first equation establishes a self-consistent relation 
between the density of q-particles p q and their chemi¬ 
cal potential p q . The above densities can be written as 
regular Fermi integrals by shifting the chemical potential 
according to p! q = p q — U q . The Fermi-Dirac distribution 
indeed reads : 

Mp) 1 + exp(^(p 2 /2m* - ^)) 


= 2 


ln(l 


- 0(- 2 


) — dp = 


h 2 

3m* 




( 21 ) 

At the thermodynamic limit, the system volume V di¬ 
verges together with the particle numbers (N q ), and the 
thermodynamics is completely defined as a function of 
the two particle densities p n ,Pp- 

We can now use the maximum-entropy variational 
principle to evaluate the mean-field approximation to the 
grand-canonical partition sum Zgc- We recall that the 
exact grand-canonical ensemble corresponds to the max¬ 
imum of the constrained Shannon entropy which is noth¬ 
ing but In Zqc '■ 


In Zgc = Sgc - (3{{H)gc - Pn{N n )GC ~ p p {N p )gc)- 

( 22 ) 

The variational principle thus states that the mean- 
field constrained entropy is the best approximation 
within the ensemble of independent-particle trial states 
to the exact maximum In Zgc'- 


In Z G c — S 0 - p((H) 0 - p n (N n ) 0 - p p (N p ) 0 ), (23) 


where the mean-field energy and particle numbers are 
defined as functions of densities, i.e. single-particle occu¬ 
pations eq. <X9, by: 

(H )o = VH ; (N q )o = Vp q (24) 


The mean-field entropy is given by : 

S 0 = \nZ 0 + /3{{W) 0 - p n {N)o - p p {N p )o) (25) 
where (W)o is the average single-particle energy : 
f°° AtTP 2 

{W) 0 = 2 J n q(p) e q(p)—f^r d P = ~ d P ln Z 0 

(26) 

with e q (p) = p 2 /2m* + U q . 

The grand-canonical partition sum in the mean-field 
approximation is thus modified with respect to the 
independent-particle partition sum Zq as : 

In Z GC * In Z 0 + (3 (<W) 0 - <£>„) (27) 
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which allows to express the grand-canonical potential 
density as a function of densities : 

9 = f5y C ~ 3 ^ + + ( CT + 1)^3 + 3 ^ e // ( 2 ®) 

At the thermodynamic limit, this quantity is equiva¬ 
lent to the system pressure P = —g. Because of ensemble 
equivalence we can then evaluate all the thermodynamic 
potentials. For example, the canonical partition sum, or 
equivalently the free energy per unit volume f, is defined 
through the Legendre transform : 

In Z 

f = - -py - = 9 + PnPn + PpPp ( 29 ) 


jjl (MeV) T = 6 MeV 

0 3 so lOO^O 



100 -50 

M Me V) 


FIG. 2: Pressure as a function of the isoscalar and isovector 
chemical potentials for a uniform system at T = 6MeV. 


FigureElillustrates the typical behavior of the pressure 
as a function of the isoscalar and isovector chemical po¬ 
tentials g = g n + g p and gz = Pn — Pp- This figure is 
computed for a uniform system at a fixed temperature. 
For symmetry reasons only the positive pz are shown. 
We can see that for some values of (p, pz), there are three 
solutions corresponding to different values for the conju¬ 
gated observables (p, p 3 ). This is the phase-transition re¬ 
gion. The true equilibrium is the solution minimizing the 
grand potential i.e. maximizing the pressure. Thus only 
the upper part of the pressure manifold corresponds to 
a thermodynamic equilibrium. At the resulting fold, the 
slope changes discontinuously, i.e. using p q = dg/dg q , 
the equilibrium uniform system jumps from a low to a 
high density solution. It is a liquid-gas first-order phase 
transition. The fold of the grand potential is more clearly 
shown in figure 0 which presents the pressure as a func¬ 
tion of g p for different values of g n . 



FIG. 3: Pressure as a function of g p for a uniform system at 
a fixed temperature T = lOMeF and regularly spaced values 
of g n from -80 to 10 MeV. 


III. PHASE COEXISTENCE 


A. Phase-coexistence conditions in multi-fluid 
systems 

For systems at the thermodynamic limit, the existence 
and the order of a phase transition are intrinsically re¬ 
lated to the singularities of the thermodynamic potential, 
— TlnZ (A) where Z (A) is the partition sum for a given 
macroscopic state characterized by the L intensive pa¬ 
rameters A = {Af} Hi- These intensive parameters are 
the Lagrange multipliers introduced in the maximization 
of the Shannon entropy under the L constraints < At > 
associated with all the L relevant observables Ai. 

A system presents a first-order phase transition if one 
of the first partial-derivatives of In Z (A) shows a discon¬ 
tinuity pj. If a non-analyticity (discontinuity or diver¬ 
gence) is present at a higher order derivative, there is 
a continuous transition. Since the equations of states 
(EO S) are the L relations < Ag > (A) = — d\ t In Z { A) 
Fl4l . we can see that a first-order phase transition corre¬ 
sponds to a discontinuity in at least one of the EOS, i.e. 
to a jump in the mean value of at least one observable at 
the transition point. This observable can then be iden¬ 
tified with an order parameter, since its value allows to 
distinguish the two coexisting phases. 

At the thermodynamic limit, the entropy S’, as all the 
extensive variables {A*,}, is additive and thus scales like 
the volume V : 

S(V, {A k }) = Vs( Pk = {A k /V}) (30) 

The volume can thus be eliminated from the thermody¬ 
namic description while the other extensive variables can 
be expressed as densities p = {p k }. If we consider two 
isolated systems of volurn Vi = aV and V 2 = (1 — a)V, 
their constrained entropy is simply the sum of the two en- 
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tropies S = S± + S 2 ■ When 1 and 2 are put into contact 
to form a system with p = ap\ + (1 — a)p 2) equilibrium is 
reached by maximization of the global entropy S so that 

S>S 1 +S 2 . (31) 

This imposes the convexity of the entropy : 

s(ap! + (1 - a)p 2 ) > as(pi) + (1 - a)s(p 2 ) (32) 

As a result, if the homogeneous system has a con¬ 
strained entropy with a convex region, a linear inter¬ 
polation between two densities pa and pb corresponds 
to a physical phase mixing of the two associated states. 
Such linear interpolations lead to the construction of 
a concave envelope which maximizes the entropy func¬ 
tional and suppresses the curvature anomaly (Gibbs con¬ 
struction). The straight lines corresponding to phase 
coexistence are defined by two points of same tangent 
plane, i.e. with identical first derivatives or intensive 
variables \ k = d Pk s , and equal values of the constrained 
entropy, i.e. equal distances between the entropy s (p) 
and the plane iZk^kPk = 0. Using [dvS{V, {A k })\/V = 
s ({/°fc}) — Sfc ^kpk this latter condition can be interpreted 
as the equality of the two system pressures. This equality 
of all intensive variables defines the conditions of phase 
equilibrium. 

B. Gibbs and Maxwell constructions 

For nuclear matter at a given temperature, the A k 
are the proton and neutron numbers: N = p n V and 
Z = p p V. Finding two points in equilibrium means find¬ 
ing two sets of densities {Pn,P p } and {Pn,P p } which 
fulfill the 3 equations Mn = Mn, P p = P p , P A = P B > 
the equality of the temperatures being insured by the 
use of an isothermal ensemble. The path corresponding 
to a constant p q (e.g. p n ) is a curve in the (p pi P) plane 
that can be determined numerically. The problem is now 
reduced to two dimensions: if this path shows a cross¬ 
ing point, there are two sets of extensive observables for 
which the intensive parameters are all equal, which is 
the condition for coexistence. This is illustrated in the 
central part of figure 0] 

From the thermostatistics point of view, using the set 
of state variables (/3,p, n ,p p ) corresponds to defining a 
neutron-grand-canonical but proton-canonical ensemble, 
denoted as pC in the following. The associated potential 
per unit volume, G p c/V = g p c, is given by : 

9pC(P, Pn: Pp) = 'M. Pnpn Pni Pp) /P (33) 

where s = S/V is the entropy density. g p c is linked by 
Legendre transforms both to the grand-canonical poten¬ 
tial gcc = —P (pressure) and to the canonical potential 
/ (free energy) : 

9pC = P T ppPp = f PnPn ■ (34) 


P (MeV.fm 3 ) G pC /Z (MeV) 


-0.1 0 0.1 0.2 0.3 0.4 -35 -30 -25 -20 -15 -10 



FIG. 4: Illustration of a Maxwell-Gibbs construction in the 
proton-canonical and neutron-grand-canonical ensemble at a 
fixed neutron chemical potential p n = —30 MeV and temper¬ 
ature T = 10 MeV . The bottom part presents the thermody¬ 
namic potential density G p c /V as a function of the proton 
density p p . The straight line shows the convex envelop inter¬ 
polating between the two phases A and B. This corresponds 
to the Maxwell construction on p p as a function of p p (left 
center). The top-right part shows the convex envelop of the 
thermodynamic potential per particle G p c/Z as a function 
of the inverse of the proton density 1 /p p which is associated 
with a Maxwell construction for the pressure as a function of 
1 /p p (top center). The two Maxwell constructions correspond 
to the crossing of the P-p p diagram (central figure). 

It is important to stress that ensembles are equivalent 
at the thermodynamic limit which is used in the pre¬ 
sented mean-field approach. The difference comes from 
the thermodynamical variables and potentials that are 
used. Therefore, the construction of the phase transition 
is not the same. The proposed ensemble leads to a sim¬ 
ple 1-dimensional Maxwell construction which is much 
simpler than the usual multidimentional Gibbs construc¬ 
tion in a multi-fluid system. Namely, for the considered 
system: 

• In the canonical ensemble, the variables are the 
densities ( p n and p p ) and the potential is the free 
energy f. The Gibbs construction is then a complex 
2 -dimensional construction of its convex envelope. 

• In the grand-canonical ensemble, the variables are 
the chemical potentials (p n and p p ) and the poten¬ 
tial is the grand potential g. In a first-order phase- 
transition region this function is multivalued. The 
Gibbs construction is then a complex 2-dimensional 
construction of the intersection of two of the three 
manifolds. 

• In the proposed ensemble, the variables are one 
density (e.g. p p ) and one chemical potential (p, n ) 
and the potential g p c is the corresponding Legen¬ 
dre transform of / and P. Since this ensemble has 
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only one density left, the phase coexistence is re¬ 
duced to a simple problem like in a 1-component 
fluid. At constant p n , g p c as a function of p p should 
be replaced by its convex envelope. This is now a 
1 -dimensional problem which is in fact the usual 
Maxwell construction. 

Numerically, we can directly perform this Maxwell con¬ 
struction on the function p p (p p ) for constant (3 and p n . 
The comparison between the results obtained using this 
method and the method described in the previous sec¬ 
tion using the p p — P crossing point gives an estimation 
of our numerical error, which comes out to be less than 
0 .2% for all temperatures. 

We can additionally remark that the above reasoning 
also holds if we study the pC-potential per proton leading 
to a Maxwell construction for the pressure P as a func¬ 
tion 1/p p . Both constructions are illustrated on figure 0] 
The convex envelop of G p c/V ( G p c/Z ) corresponds to 
an equal-area Maxwell construction on p p (p p ) (P(l/p p )) 
and to the crossing point between the two phases in the 
p p versus P graph. 

It should be noticed that the introduction of a statis¬ 
tical ensemble in which only one density is kept, all the 
other ones being replaced by their associated intensive 
parameters, can always be applied to study a phase tran¬ 
sition with a one-dimensional order parameter. The only 
condition is that the considered density is not orthogonal 
to the order parameter, i.e. its value is different in the 
two phases. In such an ensemble, the multidimensional 
Gibbs construction reduces to a simple one-dimensional 
Maxwell construction. 


IV. RESULTS 
A. Coexistence region 

In order to explore the phase diagram, the method il¬ 
lustrated in the previous section has to be performed for 
various temperatures and chemical potentials. Figure 0 
illustrates the p n dependence for T = 10 MeV. We can 
see that for a broad range of p n the p p {p p ) equations of 
state present a back-bending associated with an instabil¬ 
ity which must be corrected using a Maxwell construc¬ 
tion. This defines the transition points Pp(p n ,T) and 
P t (p n ,T). The ensemble of the transition points consti¬ 
tutes a coexistence curve at the considered temperature. 
This curve is limited by two critical chemical potentials, 
Pn(T) and p> ( T ), which in turn define the proton critical 
chemical potentials p p (T) = Pp(p<(T),T) and p p (T) = 
p p (Pn (T), T), since p p is maximum when p n is minimum. 
The transition is observed only in a finite range of tem¬ 
peratures below a given temperature T c which is nothing 
but the critical temperature of symmetric matter. 

The Gibbs construction of phase coexistence leads 
to well defined partition sums fulfilling the thermody¬ 
namic stability requirement. The resulting pressure P = 



P p (fm' 3 ) 


FIG. 5: Maxwell construction in the p p versus p p diagram for 
a fixed temperature T = 10 MeV and regularly spaced values 
of p n from —80 to 10 MeV. 

T\hZgc(T, p n , p p )/V at the temperature T = 10 MeV 
is shown in figure El 



FIG. 6: Equilibrium pressure computed at a temperature 
T = lOMeV as a function of p n and p p after performing a 
Gibbs construction for different values of p n ■ The resulting 
fold line is the coexistence line ending at two critical points. 
The traced path corresponds to a transformation at constant 
proton fraction Z/A = 0.3 (see Section IV-B). 

Along the coexistence line p p (p n ,T), the pressure 
presents a fold, i.e. the derivative perpendicular to the 
line is discontinuous. It is by definition a region of first- 
order phase transition. The coexistence line at fixed tem¬ 
perature is limited by two points {p^{T),p p {T)) and 
(p> (T), Pp(T)). They correspond to vertical tangents in 
the grand-potential first derivatives p q (p n , p P ), which are 
singularities in its second derivatives. Hence, the limit¬ 
ing points are second-order critical points, i.e points of 
continuous transition. 

We have represented in figure 0 the first derivative 
of In Zgg in the p n direction, i.e. the neutron density 
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FIG. 7: First derivative of the grand-canonical partition sum, 
pn = dp/dpn, as a function of p n and p v at a temperature 
T = lOMeV. The path Z/A = 0.3 is also shown. 


Pn = d Pn p, which is discontinuous on the first-order line 
and continuous with a vertical tangent at the two critical 
points. Because of the exact isospin symmetry of the SLy 
interaction, the proton density p p is symmetric to p n : it 
is the same surface with inversion of the axes p n and p p . 

The first derivatives of In Zqc on both sides of the 
phase transition line define pairs of points (p n , p p ) in co¬ 
existence respectively at low (gas) and high (liquid) den¬ 
sity. These phases merge together at the critical points. 
The coexistence region in the proton and neutron-density 
space is shown on figure |H] for a fixed temperature. The 
solid lines in this figure give several iso-p n paths. Be¬ 
cause the construction of the concave envelop of the con¬ 
strained entropy is nothing but a linear interpolation be¬ 
tween two phases, inside the coexistence region the iso- 
p n lines are straight lines in (p n , p p ) representation. The 
bottom part of the figure shows the coexistence region 
in the total-density and proton-fraction plane ( p,Z/A ), 
which consists in the change of variables p = p n + p p , 
y = Z/A = p P /(p n + Pp)- Because of the non-linearity 
of the variable change, we can see that coexistence does 
not correspond to a straight line in this representation, 
and the value of Z/A evolves when passing from the 
dense phase to the diluted phase in coexistence. The 
only exception is the case of symmetric nuclear matter, 
at y = Z/A = 0.5. 

The proton-fraction difference in the two phases can 
be appreciated from figure [5] which shows Z/A as a func¬ 
tion of p n and p p for a fixed temperature T = lOMeV. 
Z/A being a combination of the two order parameters p n 
and p p , it also presents a discontinuity at the first-order 
phase-transition border, with the only exception of the 
symmetric nuclear matter point where Z/A = cst = 0.5. 
Correlating this plot with figure Q one can see that the 
dense phase (e.g. point B’ or B) is systematically closer 
to isospin symmetry Z/A = 0.5 than the diluted one (e.g. 
point A or A 1 ). This phenomenon is known as isospin dis¬ 
tillation (J, 20U21]. 




FIG. 8: Closed curves : coexistence border in the pro¬ 
ton and neutron-density plane (upper part) and total-density 
and proton-fraction plane (lower part). Black dots : critical 
points. Full lines : constant p n paths for regularly spaced val¬ 
ues between p n = 10 MeV and —75 MeV. Paths AA' and BB' 
refer to a transformation at Z/A = 0.3 (see Section IV-B). 


B. Transformation at constant Z/A 


In low-energy heavy-ion collisions, the proton and neu¬ 
tron numbers obey two independent conservation laws, 
implying that the proton fraction Z/A is conserved in the 
reaction. It is therefore of interest to consider a transfor¬ 
mation at constant Z/A nm In the ( p n ,p p ) plane (or 
equivalently in the ( ptot,Z/A ) representation) Z/A = cst 
transformations are straight lines which cross the differ¬ 
ent constant-p,„ curves. Inside the coexistence region, the 
system with a given proton fraction is decomposed into 
two phases, located at the intersections of the coexistence 
































FIG. 9: Z/A = p p /{pn + pp) as a function of p n and p p for 
a fixed temperature T = lOMeV. The path at Z/A = 0.3 is 
also shown. 

curve with the corresponding constant-/r n curve. Since 
the constant-/z n curves are not aligned on constant-Z/A 
lines (except for symmetric matter), the constant-Z/A 
transformation does not make a transition from liquid 
to gas at a unique value of Hn but shows a continuous 
smooth evolution of the intensive parameters along the 
coexistence line. 
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FIG. 10: Projection of the coexistence region in the Z/A 
versus p n plane. Black dots are critical points, and stars are 
points of maximal asymmetry. The path at Z/A = 0.3 is also 
shown. 

The transformation Z/A = 0.3 is given as an example 
by the dotted lines in figure^ and the grey path in figures 
MM Let us follow this transformation from the low 
density phase. When the system reaches the coexistence 
border (point A), a liquid phase appears in B' at the same 
value of p , n , /i p and T. This point coincides with A in the 
representation of figure El We can see from figures El and 
El that the liquid fraction is closer to symmetric nuclear 
matter than the original system, as expected from the 
isospin-distillation phenomenon. 

Following the system inside coexistence along the line 
A — B, the neutron chemical potential increases. Phase 
separation implies that, in this region, the system is com¬ 


posed of two phases located at the coexistence boundaries 
and correponding to the same value of the intensive pa¬ 
rameters. The diluted phase goes along coexistence from 
A to A’ while the dense phase goes on the other side of 
the coexistence border from B’ to B. When it reaches B 
the gas is entirely transformed into a liquid, the phase 
transition is over and the Z/A = 0.3 transformation cor¬ 
responds to a homogenous system again. 

The evolution of p n during the transition can be quan¬ 
titatively discussed on figure [TUI which gives a projection 
of figure El on the fi n axis. This figure clearly shows that 
a constant-Z/A transformation does not cross the coex¬ 
istence at a unique p n value but explores a finite range of 
chemical potentials. The system is thus forced to follow 
the coexistence line in the intensive parameter space, as 
shown in figure El and in figure fTTI The only exceptions 
are symmetric matter and the 2 maximum asymmetries 
of the coexistence region for each temperature. 



FIG. 11: Coexistence line, i.e. line of first-order phase tran¬ 
sition, in the intensive-parameter plane (/i n ,/tp). The path at 
Z/A = 0.3 is also shown. 

It is important to notice that the behavior shown in 
figures m and HD is the generic behavior expected when 
a conservation law is imposed on an order parameter 
I 22 I [l23j. Indeed, the usual discontinuity of the order 
parameter characteristic of a first-order transition is pre¬ 
vented by the constraint. If the system reaches coexis¬ 
tence, the only way to fulfill the conservation law on the 
order parameter is to follow the coexistence line until 
the conservation law becomes compatible with a homo¬ 
geneous phase. 

This continuous evolution by phase mixing hides the 
EOS discontinuity associated with the first-order phase 
transition.This is illustrated by the evolution of p n and P 
as a function of the total density ptot at Z/A = 0.3 in fig¬ 
ure El Full lines are obtained after the thermodynamic 
potential is made convex by phase mixing, i.e. with the 
Gibbs construction. These functions present no plateau 
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FIG. 12: Transformation at constant Z/A = 0.3 in nuclear 
matter at T — 10 MeV. Top : pn(ptot). Bottom : P(ptot). 
Dotted lines : solution for an homogeneous system. Full lines 
: Gibbs construction. Points A and B give the coexistence- 
region borders. 

in the coexistence region, and the transformation looks 
like a continuous transition. Yet by definition the system 
is going through a first-order phase transition since the 
first derivatives of the grand potential are discontinuous. 

This clearly stresses the fact that the behavior of spe¬ 
cific transformations should not be confused with the in¬ 
trinsic thermodynamic properties. Indeed a transforma¬ 
tion is an explicit or implicit relation between the thermo¬ 
dynamic variables, e.g. a constant proton fraction Z/A 
implies p p = p n Z/(A — Z). Then the variation of the 
intensive variables conjugated to one of the order param¬ 
eters, e.g p n for p n is given by : 

(dp n _ 9 p n ^ dp n Z (35) 

V dpn J Z/A=cst ®Pn 9 Pp A — Z 

since 


{Pn{Pn))z/A=cst ~ Qp yPn,Pn^ _z) 

Then a first-order phase transition characterized by 
dp n /dp n = 0 leads to a continuous evolution of p n as a 
function of p n according to eq. (EI51) because of the con¬ 
straint of constant proton fraction. 

The isospin degree of freedom does not change the or¬ 
der of the nuclear liqu id-g as phase transition as claimed 
in different articles 0, 1131 ■ It remains first order. Only 


the constant-proton-fraction transformations (or other 
transformations constraining an order parameter) mimic 
a continuous transition because they do not cross the co¬ 
existence line at a single point, but are forced to follow 
it to fulfill the conservation law. 

In general, transformations involving a constraint on 
an order parameter always appear continuous even in the 
presence of a first order phase transition |22| H3 ■ 


C. Phase diagram: T = 0 specificity 

Until now we have presented a study at a fixed finite 
temperature. We will now consider the effect of temper¬ 
ature on the phase diagram, from the particular case of 
zero temperature to the symmetric-matter critical tem¬ 
perature T c above which there is no transition any more. 

The specificity of the zero-temperature case is the pos¬ 
sibility to reach a vanishing density with a finite chemi¬ 
cal potential l2~il| , while at any finite temperature a given 
density can be zero only if the associated chemical po¬ 
tential goes to — oo (see eq. O)- This is a trivial con¬ 
sequence of the specificity of Fermi-Dirac distribution at 
T = 0. In this case, the associated thermodynamic po¬ 
tential presents at zero density an edge with a finite slope, 
the associated finite chemical potential. An example is 
given by the (free) energy as a function of total density for 
symmetric nuclear matter in figure ^ Then, if the ther¬ 
modynamic potential for the uniform system presents a 
concave region reaching zero density, the construction of 
the convex envelope does not reduce to the usual tangent 
construction between two points in coexistence. Indeed 
the interpolating plane defined by phase mixing will not 
be tangent to the thermodynamic potential on the zero- 
density edge. 

The coexistence region at T = 0 is shown in figure 
El It can be divided into three zones corresponding to 
three kinds of equilibria with different conditions on the 
intensive parameters. In two small regions, labeled 2 in 
figure mi and associated by isospin symmetry, both pro¬ 
ton and neutron densities are finite in the gas phase. The 
usual bi-tangential construction can then be performed. 
This imposes the standard equality between all intensive 
parameters in the two phases A and B: /3(A) = (3(B), 
Pn(A) = p n (B), p p (A) = p p (B) and P(A) = P(B). 
This region is limited on one side by a critical point be¬ 
yond which there is no more curvature anomaly, and on 
the other side by the vanishing of one of the two densities 
in the gas phase. This second case corresponds to two 
symmetric regions denoted 1 in figure El There, we have 
equilibria for which the low density phase is on an edge 
p q = 0, the second density being finite. For simplicity, 
let us take the case of globally neutron-rich matter where 
p p = 0 and p n > 0 at the low-density coexistence edge. In 
this case, the convex envelop is tangent to the thermody¬ 
namic potential of the uniform system only in the dense 
phase. This is a mono-tangential construction. It means 
that phase equilibrium does not require the equality of 
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FIG. 14: Specificity of T = 0 : Mono-tangential construction 
for a fixed neutron chemical potential p n = 5 MeV. Upper 
part : corresponding path in the density plane. It reaches the 
edge p p = 0 (with a finite p „). Middle part: thermodynamic 
potential density g v c (see section III-B). The concave region 
is corrected by phase mixing, which corresponds here to a 
mono-tangential construction (straight line). Lower part: p p . 
The mono-tangential construction on g v c does not correspond 
to the usual Maxwell construction, that can not be performed 
because of the finite value of p p on the edge. 


FIG. 13: Coexistence region at T=0. It is split into 3 differ¬ 
ent types of equilibria between liquid and gas (see text). Top 
: coexistence border (thick line) in the density plane. The 
straight lines relate two coexisting phases. Bottom : coexis¬ 
tence line in the plane of chemical potentials. 


the potential derivative in the p v direction, i.e. the equal¬ 
ity of the proton chemical potentials in the two coexisting 
phases. Since the anomaly involves only the derivative in 
the p p direction, the other equilibrium conditions charac¬ 
terizing the Gibbs construction are still satisfied in this 
region, namely the neutron chemical potential, pressure 
and temperature have to be the same in the two phases. 

Since the construction involves two points with the 
same p, n , we can use the neutron-grand-canonical proton- 
canonical potential G p c (see section III-B) in order to re¬ 
duce the problem to a one-dimensional case as illustrated 
in figure ITU This is a way to solve the problem in prac¬ 


tice. Looking at a fixed p, n for which the proton density 
reaches its edge p p = 0 at zero temperature (upper part), 
we can see that a standard equal-area Maxwell construc¬ 
tion is not possible in this case (lower part). However, the 
concave region in the thermodynamic potential (middle 
part) has to be corrected by phase mixing. Equilibrium 
is then given by a mono-tangential construction on g p c- 

Such equilibrium between a neutron gas and a two- 
fluid liquid leads to a discontinuous change of p p in the 
gas. Because of the dominance of region 1 with respect 
to region 2 in the phase diagram, the simplification is 
often made in the literature |24l that phase coexistence in 
neutron-star matter can be modelized as the equilibrium 
between neutron-rich nuclei and a pure neutron gas. It 
should however be noticed that this is possible only if the 
temperature is exactly zero, a finite proton fraction being 
associated to the gas phase at any finite temperature. 

The last case, denoted 0 in figure ITT! corresponds to 
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both p p = 0 and p n = 0 in the low-density phase. This 
case corresponds to a dense phase in equilibrium with 
the vacuum, i.e. at zero pressure. This region has a 
very simple physical interpretation. If the gas phase is 
given by the point (p r , p p ) = (0,0), this means that the 
coexistence lines of zone 0 are constant Z/A lines. The 
coexistence border on the liquid side is the locus of zero 
pressure in the Z/A interval corresponding to region 0. 
For each value of Z/A, the liquid border is then given 
by the minimum of the energy per particle computed 
for this constant Z/A. Zone 0 corresponds then to the 
chemical potential interval in which a self-bound liquid 
can be defined. 


D. Phase diagram: temperature dependence 

In order to determine the phase diagram of nuclear 
matter, we establish the ensemble of points at equilib¬ 
rium in the (p n ,p p ) plane for different fixed values of 
temperature. The result is represented in figure as 
a function of the variables (p, Z/A) in order to under¬ 
line the role of isospin. Since protons and neutrons play 
symmetric roles, the resulting curves are symmetric with 
respect to the axis Z/A = 0.5. For this reason, in the 
following, only the neutron-rich part [Z/A < 0.5) will be 
discussed. 



FIG. 15: Coexistence region in the (p, Z/A) plane at different 
temperatures. Thick line : T = 0. Dotted lines : T = 4, 6, 
8MeV. Solid lines : T = 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 
1AMeV. Black dots : critical points at each temperature, 
including T c = lA.5AMeV. Dashed-dotted line : points of 
maximum asymmetry for coexistence at each temperature. 

The widest coexistence region corresponds to zero tem¬ 
perature. Equilibrium points that belong to the edge 
p p = 0 make a line at Z/A = 0. When temperature is 
introduced, coexistence region concerns non-zero values 
of p and Z/A , which means that equilibrium is always 
between phases containing both kinds of particles. 

The critical points are also reported for each temper¬ 
ature. They are second-order-transition points. They 
form a critical line. An interesting feature is that the 



FIG. 16: Cuts of the coexistence region as a function of the 
total density for regularly spaced values of Z/A between 0.05 
and 0.5. Thick grey line : line of critical points. Dashed- 
dotted line : line of maximum asymmetry. 



FIG. 17: Critical temperatures as a function of the system 
Z/A (thick grey line). The curve of maximal temperatures 
for coexistence is also shown (dashed-dotted line). 


critical density increases with the asymmetry while the 
critical temperature decreases, as can also be seen from 
figures EH and rm This is in agreement with previous 
studies using different effective interactions 0, J21| . 

It should be stressed that, on each coexistence curve 
at a fixed temperature, the critical points do not corre¬ 
spond to the minimum (or maximum) Z/A. This point 
((Z/A) mm , p mm ) is situated at a density lower than the 
critical point, i.e. on the gas side of coexistence. As a re¬ 
sult, transformations can be performed at constant Z/A 
such that the system enters at a point of coexistence as 
a gas and exits at another point still as a gas. This is 
the so-called retrograde transition [l|]. This happens for 
the small region for which ( Z/A) mm < Z/A < Z/A c . 
Between those two points, there is appearance and dis¬ 
appearance of a liquid phase at ( Z/A) L > Z/A at equi¬ 
librium with a gas at ( Z/A) G < Z/A. For a too neutron- 
rich system such that Z/A < (. Z/A ) mm , there can be no 
phase coexistence. 
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As temperature grows, the coexistence region is re¬ 
duced. The values of (Z/A) c and ( Z/A) mm become closer 
to the symmetry Z/A = 0.5. Densities p c and p mm fol¬ 
low opposite evolutions: p min grows with temperature, 
while p c diminishes. The coexistence region disappears 
at the critical temperature T c , for which only symmetric 
matter presents a second-order transition point. Criti¬ 
cal and minimum Z/A lines join at this ultimate critical 
point. 

It is also interesting to look at the coexistence zone 
as a function of Z/A (see figures ITTH IT71 . We can see 
that the dependence of the critical temperature (as well 
as the maximal temperature) on Z/A is weak. Only for 
very high values of asymmetry the difference between the 
two temperatures for a fixed Z/A can be several MeV. 
This is another way to visualize the retrograde transition, 
meaning that we can have saturated vapor at supercriti¬ 
cal temperatures. 



(MeV) 

FIG. 18: Phase diagram in {pn,p P ) for different values of T. 
Thick line : T = 0. Dotted lines : T — 4, 6, 8 MeV. Solid 
lines : T = 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, UMeV. 
Black dots : critical points at each temperature, including 
T c = 14.54 MeV. 

Finally, it is interesting to look at the coexistence man¬ 
ifold in the (/3,p n ,p p ) space, as shown in figure^] One 
can see that the coexistence is almost perpendicular to 
the ps = 0 axis, stressing the fact that the nuclear liquid- 
gas phase transition is dominated by its isoscalar compo¬ 
nent @,01- 

E. Critical behavior 

Let us now study in more details the critical behav¬ 
ior of the system. Conversely to the usual single-fluid 
liquid-gas phase transition which presents a unique crit¬ 
ical point, a two-fluid system is critical along a line in 


the (/?, p n , p p ) intensive-parameter space. The criti¬ 
cal points are obtained for each temperature (below the 
symmetric-matter critical temperature) by determining 
the p n value for which the two phases with the same 
tangent plane merge together. To evaluate this point, we 
have used both the Maxwell construction in the neutron- 
grand-canonical proton-canonical ensemble (see section 
III-B) and the existence of a crossing point in the P ver¬ 
sus p p curve at constant (T, p n ). The critical points also 
correspond to the disappearance of a concave region in 
the uniform-system free energy. This disappearance of 
the spinodal region provides an independent way to eval¬ 
uate the location of the critical line: we can represent as 
a function of p n the lower value of the free-energy cur¬ 
vature in the (p ni p p ) plane. The resulting curve crosses 
the horizontal axis at the critical value of p n . Below the 
symmetric-matter critical point T c , these methods lead 
to the definition of two symmetric critical lines, one for 
neutron-rich systems (p n = p> (T), p p = p p (T) ) and the 
other for the opposite isospin (p p = p p (T), p n = p^{T) 



log(5) 


FIG. 19: Solid lines with points : evolution of p 3 — pz c as 
a function of S (see text) near the critical point for different 
temperatures : T = 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 
14 MeV. Straight dashed-dotted lines : power laws with crit¬ 
ical exponent /3p = 2. 

At the approach of the critical point, the distance in 
the space of observables between the two phases (A and 
B) in equilibrium, i.e. the order parameter, goes to zero 
as a power of the distance to the critical point in the 
intensive-variable space. The resulting power law is char¬ 
acterized by the critical exponent (3. In order to study 
this behavior for nuclear matter in the (p n , p p ) plane, we 

consider a distance 5 = yj(p& — p£) 2 + (p^ — Pp) 2 - In 
a two-fluid system the behavior of this distance can be 
studied in two ways. 

First, for a fixed temperature T = /3 _1 , 6 is expected 
to depend on the chemical potentials as : 

Sf3(p q ) oc ( p c q - /r 9 ) 1 /^ 


(37) 
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FIG. 20: Points : evolution of T — T c as a function of 8 (see 
text) near the critical point for different fi n : fi n = fif(T c ), 
Hq{T = 12 MeV), fif(T = 10 MeV). Straight lines: power 
laws with critical exponent /3 Mn = 2. 

where f3p is the critical exponent at a fixed temper¬ 
ature and fi c q the critical value of fi q for the considered 
temperature. Since fi n and fi p are constrained to be on 
the coexistence line, the above relation is independent 
of the chemical potential selected to perform the study. 
For symmetry reasons we have introduced the isoscalar 
fi = fi n + Mp and the isovector [13 = fi n — Up chemical 
potentials. Since the dependence on fi is negligeable (see 
figure EE3, we have focused our study on the evolution 
of S with fi 3 . The results are presented in figure ITTT1 We 
can observe that our results perfectly follow the expected 
critical behavior with the mean-field value (3p = 2. 

For a fixed value of fi q , the distance S Pq goes to zero as 
the temperature tends to the critical one (3 C . The power 
law should be : 

6^(f3) cc ({3 C - (38) 

where P llq is the critical exponent at a fixed fi q . Again, 
since the temperature (3 and the chemical potential fi q in 
coexistence are related by a Clapeyron-like relation, the 
above study leads to the same critical behavior if studied 
as a function of fi q > with q' ^ q, instead of (3. In figure [20 
we present the evolution with temperature for different 
chemical potentials. This graph shows that our results 
fulfill the mean-field scaling (3 tlq = 2. 

V. CONCLUSION 

In this paper, based on a mean-field analysis of nuclear 
matter with a realistic Skyrme SLy230a effective interac¬ 


tion, we have established that nuclear matter presents 
a first-order phase transition even when the isospin de¬ 
gree of freedom is explicitly accounted for. This results 
from the existence of a spinodal region, which is a re¬ 
gion where the free energy of a homogeneous system is 
concave. In the case of infinite systems, such curvature 
anomaly is corrected by constructing the convex enve¬ 
lope. This is a tangent construction that links pairs of 
points in the space of observables with the same values 
for the intensive parameters. It corresponds to points 
of discontinuous first derivatives for the grand-canonical 
potential of the system in the space of Lagrange inten¬ 
sive parameters, along a coexistence manifold limited by 
a critical line. Except on this limit which corresponds to 
a second-order phase transition, the slope discontinuity 
demonstrates that the system is undergoing a first-order 
phase transition. For fixed values of temperature (below 
the symmetric-matter critical one), coexistence lines are 
obtained in the chemical-potential plane. They are lim¬ 
ited by critical points that correspond to proton fractions 
depending on the temperature. As lower temperatures 
are considered, more asymmetric nuclear matter can be 
involved in a first-order phase transition. 


As for the study of critical behaviors, we have found 
that all the numerical data can be fitted by power laws 
with critical exponents equal to 2 ? which is consistent 
with generic mean field predictions j22j . Calculations be¬ 
yond mean field are needed to obtain a result which could 
be characteristic of a given universality class. Looking at 
the isospin content of the phases, we show that the proton 
fraction is discontinuous at the phase transition except 
at the critical points and for symmetric matter. In asym¬ 
metric nuclear matter, the proton fraction can be used 
as an order parameter. When the transition occurs, the 
liquid gets closer to symmetry while the gas is enriched 
in the more abundant species. This is the well-known 
isospin fractionation. This has a strong influence on the 
constant-proton-fraction transformations, since in order 
to fulfill the imposed conservation on an order parame¬ 
ter, the transformation is forced to follow the coexistence 
line instead of crossing it. This hides the slope disconti¬ 
nuity characteristic of a first-order phase transition, the 
transformation appearing as continuous. 
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